% Process the eej strength derived from satellite data
% and compare it with ground deribed eej indices
%latest date 30 OCT 2005
clear;
data = load('c:\manoj\projects\eej\eej_par_all.txt');
%
%data1 = load('c:\manoj\projects\eej\eejdata.txt'); %<ettmnhyb 2000-2001; 
%data1 = load('c:\manoj\projects\eej\eejdata_new.txt');%<ettmnhyb2000-2002
%and corrected tirmnabg

%----31 Jan 2006 Added HUA - FUQ test 
%data1 = load('c:\manoj\projects\eej\eejdata31jan.txt');
%mjd,y,m,d,j-1,ettmnhyb(pp),ettmnabg(pp),tirmnhyb(pp), tirmnabg(pp),huamnfuq(pp),kp);

%-------------------------------------------------------

%---1Feb added MBO-GUI
%data1 = load('c:\manoj\temp\eejdata_all.txt');
%data1 = load('c:\manoj\temp\eejdata_single.txt');
%15 Jan 2006 Some notes
% Basically tried to veryfy the results. Checked the time stamps ->
% corrected the IST-Solar Time problems.. Checked with many LT bin
% combinations
% someof the resluts are mentioned in correlation.doc 
% LT 9-15 and bins 10 degree with the necessary corrections brings the most
% aestetic (?) results


%-------- correlate the CHAMP derived EEJ strength and ground delta H
% with LT correction
%--------------------------------------------------------------------
bin5 = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
bin10 = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
bin20 = [-60 -40 -20 0 20 40 60];

%Control Center

obsdata = 9;  %ETT-HYB = 1,ETT-ABG = 2;TIR-HYB = 3;TIR-ABG = 4; HUA-FUQ 5; MBO-GUI = 6; AAE-QSB 7; GUA-CBI 8; YAP-BIK 9;
satdata = 2; % %Pear Current Density = 1; av cur distan 2;
obslt = [9.9 13.1]; %NOTE the LT are represented as for eg 8.99999 or 9.001 the filter
incrlat = 10;
x1data = bin10;
mindate = tmjd(2000,01,01,0,0,0,1); % starting date
maxdate = tmjd(2002,12,31,0,0,0,1)+1; % starting date
sqcor = 0;
plotopt = 'g.-';
plotcol = 'g';

%  index = [];
%  save c:\manoj\temp\index index;
%load c:\manoj\projects\eej\guacm4.mat index sq_est

switch obsdata,
    case 1,
        load c:\manoj\projects\eej\EEJSingle ETTMNHYB;
        k = nanmean(ETTMNHYB);
        k=circshift(k,[0,5])'; % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = data1(:,1) + ((78-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'b*-';

    case 2,
       load c:\manoj\projects\eej\EEJSingle ETTMNABG;
        k = nanmean(ETTMNABG');
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 8;
        data1(:,1) = data1(:,1) + ((78-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
        obslon = 82.5;
         st = 12.8;
        en = 132.8;
        plotoptions = 'b*-';

    case 3,
       load c:\manoj\projects\eej\EEJSingle TIRMNHYB;
        k = nanmean(TIRMNHYB');
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 9;
        data1(:,1) = data1(:,1) + ((77.85-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78.5 (TIR) and 82.5 (Allahabad - IST)
        xtime = [1:24]';
         obslon = 82.5;
         st = 12.8;
        en = 132.8;
        plotoptions = 'b*-';

    case 4,
        load c:\manoj\projects\eej\EEJSingle TIRMNABG;
        k = nanmean(TIRMNABG);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 10; 
        data1(:,1) = data1(:,1) + ((75-82.5)/15)/24; % This assumes 75 (between TIR and ABG)
        xtime = [1:24]';
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'k.-';

    case 5,
        load c:\manoj\projects\eej\EEJSingle HUAMNFUQ;
        k = nanmean(HUAMNFUQ,2);
        k = circshift(k,19);
        xtime = 0.5:1:23.5;       
        obs_index = 11;
        obslon = 285;
        st = -140;
        en = -10;
%         st  = -139.5
%         en = -10.5
        plotoptions = 'k*-';
    case 6,
        load c:\manoj\projects\eej\EEJSingle MBOMNGUI;
        k = nanmean(MBOMNGUI)';
        k = circshift(k,23); % 23 ~ = 343/15
        xtime = 0.5:1:23.5;       
        obs_index = 12;
        obslon = 343;
        st = -81.5;
        en = 47.5;
        plotoptions = 'g*-';
    case 7,
        load c:\manoj\projects\eej\EEJSingle AAEMNQSB;
        k = nanmean(AAEMNQSB)';
        k = circshift(k,2); % 2 ~ = 37.5/15
        xtime = 1:24;       
        obs_index = 13;   
        obslon = 37.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'k*-';
 %             plotoptions = 'r*-';
    case 8,
        load c:\manoj\projects\eej\EEJSingle GUAMNCBI;
        k = nanmean(GUAMNCBI)';
        k = circshift(k,9); % 9 ~ = 144/15 (the remaining would be taken care by shifting time axis to 1:24 (instead 0.5:1:23.5)
        xtime = 1:24;       
        obs_index = 14;  
        obslon = 143;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'r*-';
    case 9,
        load c:\manoj\projects\eej\YAPMNBIK;
        k = nanmean(YAPMNBIK);
        k = circshift(k,[1,554]); %554 = (138.5/15)*60 minutes (UT shift in minutes)
        xtime = (linspace(0,24,1440));       
        obs_index = 2;
        data1(:,1)= yapdate;
        data1(:,2) = yapmnbik;
        %obslon = 138.5;
        obslon = 138.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'c*-';
    otherwise,
        display('Error');
        break;
        return;
end;

switch satdata,
    case 1,
        eej_index = 9;
    case 2,
        eej_index = 10;
    otherwise,
        display('Error');
        break;
        return;
end;        

if sqcor == 1,
    data1(index,obs_index) = data1(index,obs_index)-(sq_est(:,1)-sq_est(:,2));
   % data1(index,obs_index) = (sq_est(:,1)-sq_est(:,2));
end;

[p,s] = polyfit(xtime,k,20);
%----------------------------

ndd = zeros([1,25]);
nloop = 1;
lto = [];
ltc = [];
xdata = [];
nbindata = [];
ceff = [];
cerr = [];
nbin = 1;

for ki = st:incrlat:en,
L = data(:,2) > ki & data(:,2) < ki+incrlat & data(:,1) < maxdate  &  data(:,1) > mindate; %The last filter is to select only data before Dec 31, 2002
IndiaE = data(L,:);
ndata = 1;
ind1 = [];
ind2 = [];
ltt1 = [];
ltt2 = [];
nnodata=1;
for i = 1:length(IndiaE(:,1)),
    ind = findnearest(IndiaE(i,1),data1(:,1));
    [lt1,dummy1] = champ_lt(IndiaE(i,1),IndiaE(i,2),obslon);
    [dummy2,lt2] = champ_lt(data1(ind,1),IndiaE(i,2),obslon); % This is because observatory FDAY ~= CHAMP FDAY
    ltt1(i) = lt1;
    ltt2(i) = lt2;
    ltt3(i) = dummy1;
    ltt4(i) = dummy2;
    obst(i) = data1(ind,1);
    chpt(i) = IndiaE(i,1);
    
         if lt2 >= obslt(1) & lt2 <= obslt(2) & abs(data1(ind,obs_index)) < 300 & data1(ind,1) < max(data(:,1)), %The last fliter take care the nodata (9999) points
            ind1(ndata) = ind;
            f_ch(ndata) = polyval(p,lt1); %Getting the expexted EEJ strength from poly coefficients for champ LT
            f_ob(ndata) = polyval(p,lt2); %Getting the expexted EEJ strength from poly coefficients for obs LT
%             f_ch(ndata) = interp1(1:24,k,lt1,'spline'); %Found spline interpolation working better 18.1.06 !
%             f_ob(ndata) = interp1(1:24,k,lt2,'spline');
%             fprintf('Lon = %5.1f i = %d f_ch = %f f_ob = %f lt1 = %f lt2 = %f kp = %3.1f\n', ki,i,f_ch(ndata),f_ob(ndata),lt1,lt2,data1(ind,end));
            ndata = ndata+1;
        else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
        end;
end;

if length(ind2) > 0,
    IndiaE(ind2,:) = [];
    obst(ind2) = [];
    chpt(ind2) = [];
end;

if length(ind1) > 1,
    obsE = data1(ind1,:);
    X = IndiaE(:,eej_index);
    Y = obsE(:,obs_index).*(f_ch./f_ob)';
    if length(X) > 3,
    [rs,rstat] = robustfit(X,Y);
    residuals = abs(rstat.resid);
    LLL = residuals == max(residuals); %Just remove the most noisy point
    X(LLL) = [];
    Y(LLL) = [];
    end;
    [cof,err] = corrcoef(X,Y);  % correlation & Correction for observatory data 

    %   [cof,err] = corrcoef(IndiaE(:,eej_index),obsE(:,obs_index));
    ceff(nloop) = cof(2,1);
    cerr(nloop) = err(2,1);
    fprintf('%4.3f ',cof(2,1));
    xdata(nloop) = ki+incrlat/2;
    ndd(nloop) = ndata-1;
    lto(nloop) = min(ltt2);
    ltc(nloop) = max(ltt2);
    nloop = nloop +1;
 % The following line is written to save the mjd indices for CM4 siumlation    
%    load c:\manoj\temp\index index;
%    index = [index ind1];
%    save c:\manoj\temp\index index;
 %-Over comment   
    %----------- PLOT SCATTER
% subplot(2,4,nloop-1);
% 
% % figure(1);
% % hold on;
% if x1data(nbin) == 0,
% plot(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)',[plotcol 's'],'MarkerFaceColor',plotcol);
% hold on;
% axis([0 0.5 0 200]);
% h=xlabel(['Longitude Bin ' sprintf('%4.1f ',ki) '^\circ  '  sprintf('%4.1f ',ki+incrlat) '^\circ']);
% set(h,'FontSize',16);
% aaa = robustfit(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)');
% [c,sfit] = polyfit(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)',1);
% m = polyval(c,[min(IndiaE(:,eej_index)),max(IndiaE(:,eej_index))]);
% %dummy = corrcoef(rstat.resid,(obst-chpt)*24);
% %plot(obst,chpt,'r.');
% % % xcor1(nloop-1) =dummy(2,1);
% plot([min(IndiaE(:,eej_index)),max(IndiaE(:,eej_index))],m,plotcol,'LineWidth', 2);
% text(0,180,sprintf('dH = %5.1f * I + %5.2f',c(1),c(2)),'FontSize',16,'Color',plotcol);
% h=title(sprintf('%d ',x1data(nbin)));
% set(h,'FontSize',16);
% end;


% 
%-----------------END
else,
    x1data(nbin) = NaN;
     
end;
nbin = nbin + 1;
clear ind f_ch f_ob ind1 ind2 ndata nnodata IndiaE obsE ltt1 ltt2 ltt3 ltt4 obst chpt;
end;

x1data(isnan(x1data)) = [];
fprintf('\n');
figure(2);
hold on;
errorbar(x1data,ceff,cerr/2,plotopt,'LineWidth',2);
axis([-80 80 -.6 1.4]);
grid on;

% %Notes
% %16.1.2006 Checked KP of all the selected events all are < 2.0 as expected
% %since the satellite data are screened for kp
% % lt2= ltt2;
% % lt2(ind2) = [];
% % L = lt2 >= 10 & lt2 <= 13;
% subplot(421);
% plot(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)','rs','MarkerFaceColor','r');
% %hold on;plot(IndiaE(L,eej_index),obsE(L,obs_index).*(f_ch(L)./f_ob(L))','bs','MarkerFaceColor','b');
% xlabel('Pear Current Density (CHAMP)');
% ylabel('\delta H nT');
% axis([0 0.5 0 200]);
% %legend(['9-15: Coef ' sprintf('%4.2f',cof(2))],['10-3: Coef ' sprintf('%4.2f',cof1(2))]);
% legend(['10-3: Coef ' sprintf('%4.2f',cof(2))]);
% % title(['ETT HYB  ' sprintf('%4.1f ',ki) '^\circ  '  sprintf('%4.1f ',ki+incrlat) '^\circ']);
% % 
% load c:\manoj\temp\ettdelta;
% L = isnan(ettdelta);
% ettdelta(L) = 9999;
% data1(:,7) = ettdelta;
